#include "constant-advection.hh"
#include <map>
#include <math.h>

namespace Dune {
  namespace Equation {
    ConstantAdvection_<1>::ConstantAdvection_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    ConstantAdvection_<2>::ConstantAdvection_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    ConstantAdvection_<3>::ConstantAdvection_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    ConstantAdvection_<1>::~ConstantAdvection_(){
    }
    ConstantAdvection_<2>::~ConstantAdvection_(){
    }
    ConstantAdvection_<3>::~ConstantAdvection_(){
    }
    void ConstantAdvection<1>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      set_equ_para("u0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
    }
    void ConstantAdvection<2>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      set_equ_para("u0",0.0);
      set_equ_para("v0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["v0"]=ptree.get("equation.v0",0.0);
    }
    void ConstantAdvection<3>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      set_equ_para("u0",0.0);
      set_equ_para("v0",0.0);
      set_equ_para("w0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["v0"]=ptree.get("equation.v0",0.0);
      para["w0"]=ptree.get("equation.w0",0.0);
    }

    double ConstantAdvection<1>::local_max_eigenvalue(double* var) {
      return std::abs(para["u0"]);
    }
    double ConstantAdvection<2>::local_max_eigenvalue(double* var) {
      return std::sqrt(std::pow(para["u0"],2)+std::pow(para["v0"],2));
    }
    double ConstantAdvection<3>::local_max_eigenvalue(double* var) {
      return std::sqrt(std::pow(para["u0"],2)+std::pow(para["v0"],2)+std::pow(para["w0"],2));  
    }

    void ConstantAdvection<1>::get_1d_fx(double* fequ, double* var){
      fequ[0] = para["u0"] * var[0];
    }
    void ConstantAdvection<2>::get_2d_fx(double* fequ, double* var){
      fequ[0] = para["u0"] * var[0];
    }
    void ConstantAdvection<2>::get_2d_fy(double* fequ, double* var){
      fequ[0] = para["v0"] * var[0];
    }
    void ConstantAdvection<3>::get_3d_fx(double* fequ, double* var){
      fequ[0] = para["u0"] * var[0];
    }
    void ConstantAdvection<3>::get_3d_fy(double* fequ, double* var){
      fequ[0] = para["v0"] * var[0];
    }
    void ConstantAdvection<3>::get_3d_fz(double* fequ, double* var){
      fequ[0] = para["w0"] * var[0];
    }

    void ConstantAdvection_<1>::get_equation_flux(int whichdim, double* fequ, double* var){
      get_1d_fx(fequ,var);
    }
    void ConstantAdvection_<2>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_2d_fx(fequ,var);
      if(whichdim==1)
        get_2d_fy(fequ,var);
    }
    void ConstantAdvection_<3>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_3d_fx(fequ,var);
      if(whichdim==1)
        get_3d_fy(fequ,var);
      if(whichdim==2)
        get_3d_fz(fequ,var);
    }
  } // Equation
} // Dune
